C************************************************************* C Test program for subroutine LIBEM3 by Stephen Kirkup * C************************************************************* C C Version 2 (Version 1 2001). Revised test problems. Changes C in LIBEM3. C School Engineering, University of Central Lancashire C www.uclan.ac.uk C smkirkup@uclan.ac.uk C http://www.researchgate.net/profile/Stephen_Kirkup C C This open source code can be found at C www.boundary-element-method.com/fortran/LIBEM3_T.FOR C C Issued under the GNU General Public License 2007, see gpl.txt C C Part of the the author's open source BEM packages. C All codes and manuals can be downloaded from C www.boundary-element-method.com C C*************************************************************** C This program is a test for the subroutine LIBEM3. The program computes C the solution to a Laplace problem interior to a sphere by the C boundary panel method. C C Background C ---------- C C We wish to solve the Laplace equation C C __ 2 C \/ {\phi} = 0 C C C For the interior problem, the domain lies interior to a closed C boundary S. The boundary condition may be Dirichlet, Robin or C Neumann. It is assumed to have the following general form C C {\alpha}(q) {\phi}(q) + {\beta}(q) v(q) = f(q) C C where {\phi}(q) is the velocity potential at the point q on S, v(q) C is the derivative of {\phi} with respect to the outward normal to S C at q and {\alpha}, {\beta} and f are real-valued functions defined C on S. C C Subroutine LIBEM3 accepts a description of the boundary of the domain C and the position of the interior points where the solution ({\phi}) C is sought, the boundary condition and returns the solution ({\phi} C and v) on S and the value of {\phi} at the interior points. C C The test problems C ----------------- C C In this test the domain is a sphere of radius 1 (metre). The solution C to the problem with a Dirichlet boundary condition ({\alpha}=1, C {\beta}=0) on the upper surface (z>=0) and with a Neumann boundary C condition ({\alpha}=0, beta=1) on the lower surface (z<0) are C sought. C Assuming the sphere is centred at the origin the boundary conditions C are specified through taking the solution to be determined by C C {\phi} = x + y + z , C C which is clearly a solution of the Laplace equation. C C The boundary is described by a set of NS=36 planar triangular panels C of approximately equal size. The boundary solution points are the C centres of the panels. C The solution is sought at the interior points (0,0), (0,0.5), C (0.25,0.25). C---------------------------------------------------------------------- C The PARAMETER statement C ----------------------- C There are four components in the PARAMETER statement. C integer MAXNS : The limit on the number of boundary panels. C integer MAXNV : The limit on the number of vertices. C integer MAXNPI : The limit on the number of interior points. C External modules related to the package C --------------------------------------- C subroutine LIBEM3: Subroutine for solving the interior Laplace C equation (file LIBEM3.FOR contains the subroutine LIBEM3) C subroutine H3LC: Returns the individual discrete Laplace integral C operators. (file H3LC.FOR contains H3LC and subordinate routines) C subroutine GLS: Solves a general linear system of equations. C (file GLS.FOR contains CGSL and subordinate routines) C file GEOM3D.FOR contains the set of relevant geometric subroutines C The program PROGRAM LIBEM3T IMPLICIT NONE C VARIABLE DECLARATION C -------------------- C PARAMETERs for storing the limits on the dimension of arrays C Limit on the number of panels INTEGER MAXNS PARAMETER (MAXNS=40) C Limit on the number of vertices (equal to the number of panels) INTEGER MAXNV PARAMETER (MAXNV=30) C Limit on the number of points interior to the boundary, where C solutions are sought INTEGER MAXNPI PARAMETER (MAXNPI=6) C Constants C Real scalars: 0, 1, 2, pi REAL*8 ZERO,ONE,THREE C Geometrical description of the boundary(ies) C Number of panels and counter INTEGER NS,IS C Number of central points (on S) and counter INTEGER NSP,ISP C Number of vetices and counter INTEGER NV,IV C Index of nodal coordinate for defining boundaries (standard unit is C metres) REAL*8 VERTEX(MAXNV,3) C The three nodes that define each panel on the boundaries INTEGER SELV(MAXNS,3) C The points interior to the boundary(ies) where the solution C is sought and the directional vectors at those points. C [Only necessary if an interior solution is sought.] C Number of interior points and counter INTEGER NPI,IPI C Coordinates of the interior points REAL*8 PINT(MAXNPI,3) C Data structures that are used to define each test problem in turn C and are input parameters to LIBEM3. C SALPHA(j) is assigned the value of {\alpha} at the centre of the C j-th panel. REAL*8 SALPHA(MAXNS) C SBETA(j) is assigned the value of {\beta} at the centre of the C j-th panel. REAL*8 SBETA(MAXNS) C SF(j) is assigned the value of f at the centre of the j-th panel. REAL*8 SF(MAXNS) C Incident field C real SIPHI(MAXNS): The incident velocity potential at the C centres of the panels REAL*8 SIPHI(MAXNS) C real PDIPHI(MAXNPI): The incident velocity potential at the chosen C interior points REAL*8 PDIPHI(MAXNPI) C Validation and control parameters for LIBEM3 C Switch for particular solution LOGICAL LSOL C Validation switch LOGICAL LVALID C The maximum absolute error in the parameters that describe the C geometry of the boundary. REAL*8 EGEOM C Output from subroutine LIBEM3 C The velocity potential (phi - the solution) at the centres of the C panels REAL*8 SPHI(MAXNS) C The normal derivative of the velocity potential at the centres of the C panels REAL*8 SVEL(MAXNS) C The velocity potential (phi - the solution) at interior points REAL*8 PDPHI(MAXNPI) C Discrete Operators (L_SS and M_SSPLUS changed in GLS) REAL*8 L_SS(MAXNS,MAXNS) REAL*8 M_SSPLUS(MAXNS,MAXNS) REAL*8 L_DS(MAXNPI,MAXNS) REAL*8 M_DS(MAXNPI,MAXNS) C Further information from system solution GLS INTEGER*4 PERM(MAXNS) LOGICAL XORY(MAXNS) REAL*8 C(MAXNS) C Work Space REAL*8 WKSPC1(MAXNS) REAL*8 WKSPC2(MAXNS) REAL*8 WKSPC3(MAXNS) REAL*8 PHI,DPHI REAL*8 SIZE3 REAL*8 P(3),NP(3),SIZEP C Counter through the x,y coordinates INTEGER ICOORD C The coordinates of the centres of the panels REAL*8 SELCNT(MAXNS,3) C INITIALISATION C -------------- C Set constants ZERO=0.0D0 ONE=1.0D0 THREE=3.0D0 C Describe the nodes and panels that make up the boundary C :The unit sphere, centred at the origin is divided C : into NS=36 uniform panels. VERTEX and SELV are defined C : anti-clockwise around the boundary so that the normal to the C : boundary is assumed to be outward C :Set up nodes C : Set NS, the number of panels NS=36 C : Set NV, the number of vertices (equal to the number of panels) NV=NS C : Set coordinates of the nodes C Set up VERTEX, the coordinates of the vertices of the panels NV=20 DATA ((VERTEX(IV,ICOORD),ICOORD=1,3),IV=1,20) * / 0.000D0, 0.000D0, 1.000D0, * 0.000D0, 0.745D0, 0.667D0, * 0.645D0, 0.372D0, 0.667D0, * 0.645D0,-0.372D0, 0.667D0, * 0.000D0,-0.745D0, 0.667D0, * -0.645D0,-0.372D0, 0.667D0, * -0.645D0, 0.372D0, 0.667D0, * 0.500D0, 0.866D0, 0.000D0, * 1.000D0, 0.000D0, 0.000D0, * 0.500D0,-0.866D0, 0.000D0, * -0.500D0,-0.866D0, 0.000D0, * -1.000D0, 0.000D0, 0.000D0, * -0.500D0, 0.866D0, 0.000D0, * 0.000D0, 0.745D0,-0.667D0, * 0.645D0, 0.372D0,-0.667D0, * 0.645D0,-0.372D0,-0.667D0, * 0.000D0,-0.745D0,-0.667D0, * -0.645D0,-0.372D0,-0.667D0, * -0.645D0, 0.372D0,-0.667D0, * 0.000D0, 0.000D0,-1.000D0/ C : Set nodal indices that describe the panels of the boundarys. C : The indices refer to the nodes in VERTEX. The order of the C : nodes in SELV dictates that the normal is outward from the C : boundary into the domain. DATA ((SELV(IS,ICOORD),ICOORD=1,3),IS=1,36) * / 1, 3, 2, 1, 4, 3, 1, 5, 4, 1, 6, 5, * 1, 7, 6, 1, 2, 7, 2, 3, 8, 3, 9, 8, * 3, 4, 9, 4,10, 9, 4, 5,10, 5,11,10, * 5, 6,11, 6,12,11, 6, 7,12, 7,13,12, * 7, 2,13, 2, 8,13, 8,15,14, 8, 9,15, * 9,16,15, 9,10,16, 10,17,16, 10,11,17, * 11,18,17, 11,12,18, 12,19,18, 12,13,19, * 13,14,19, 13, 8,14, 14,15,20, 15,16,20, * 16,17,20, 17,18,20, 18,19,20, 19,14,20/ C Set the centres of the panels, the central points DO IS=1,NS SELCNT(IS,1)=(VERTEX(SELV(IS,1),1) * +VERTEX(SELV(IS,2),1)+VERTEX(SELV(IS,3),1))/THREE SELCNT(IS,2)=(VERTEX(SELV(IS,1),2) * +VERTEX(SELV(IS,2),2)+VERTEX(SELV(IS,3),2))/THREE SELCNT(IS,3)=(VERTEX(SELV(IS,1),3) * +VERTEX(SELV(IS,2),3)+VERTEX(SELV(IS,3),3))/THREE END DO C Set the points in the domain where the solutions are sought, PINT. C : Set four points NPI=4 DATA ((PINT(IPI,ICOORD),ICOORD=1,3),IPI=1,4) * / 0.500D0, 0.000D0, 0.0000D0, * 0.000D0, 0.000D0, 0.250D0, * 0.000D0, 0.000D0, 0.500D0, * 0.000D0, 0.000D0, 0.750D0/ C The number of points on the boundary is equal to the number of C panels NSP=NS C :Switch for particular solution LSOL=.TRUE. C :Switch on the validation of LIBEM3 LVALID=.TRUE. C :Set EGEOM EGEOM=1.0D-6 C Open file for the output data OPEN(UNIT=20,FILE='LIBEM3.OUT') WRITE(20,*) 'Tests on LIBEM3' WRITE(20,*) 'A subroutine for solving the 3D interior * Laplace Equation' WRITE(20,*) 'The surface is the unit sphere that is divided * into 36 panels' C TEST 1 WRITE(20,*) WRITE(20,*) 'Test 1 - Dirichlet - {\phi}=1' WRITE(20,*) C Set up particular alpha and beta functions for Dirichlet test C and type of boundary condition DO 140 ISP=1,NSP SALPHA(ISP)=1.0D0 SBETA(ISP)=0.0D0 SF(ISP)=1.0D0 140 CONTINUE DO 160 ISP=1,NSP SIPHI(ISP)=ZERO 160 CONTINUE DO 170 IPI=1,NPI PDIPHI(IPI)=ZERO 170 CONTINUE CALL LIBEM3(MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPI,NPI,PINT, * SALPHA,SBETA,SF,SIPHI,PDIPHI, * LSOL,LVALID,EGEOM, * SPHI,SVEL,PDPHI, * L_SS,M_SSPLUS,L_DS,M_DS, * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) WRITE(20,*) 'Solutions at the interior points' WRITE(20,*) ' P(1) P(2) P(3) Computed Exact' DO 180 IPI=1,NPI P(1)=PINT(IPI,1) P(2)=PINT(IPI,2) P(3)=PINT(IPI,3) WRITE(20,999) P(1),P(2),P(3),PDPHI(IPI),1.0 180 CONTINUE C TEST 2 WRITE(20,*) WRITE(20,*) 'Test 2 - Half Dirichlet, half Neumann - * {\phi}=1, v=0' WRITE(20,*) C Set up particular alpha and beta functions for Dirichlet/Neumann test C and type of boundary condition DO 240 ISP=1,NSP/2 SALPHA(ISP)=1.0D0 SBETA(ISP)=0.0D0 SF(ISP)=1.0D0 240 CONTINUE DO 250 ISP=NSP/2+1,NSP SALPHA(ISP)=0.0D0 SBETA(ISP)=1.0D0 SF(ISP)=0.0D0 250 CONTINUE DO 260 ISP=1,NSP SIPHI(ISP)=ZERO 260 CONTINUE DO 270 IPI=1,NPI PDIPHI(IPI)=ZERO 270 CONTINUE CALL LIBEM3(MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPI,NPI,PINT, * SALPHA,SBETA,SF,SIPHI,PDIPHI, * LSOL,LVALID,EGEOM, * SPHI,SVEL,PDPHI, * L_SS,M_SSPLUS,L_DS,M_DS, * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) WRITE(20,*) 'Solutions at the interior points' WRITE(20,*) ' P(1) P(2) P(3) Computed Exact' DO 280 IPI=1,NPI P(1)=PINT(IPI,1) P(2)=PINT(IPI,2) P(3)=PINT(IPI,3) WRITE(20,999) P(1),P(2),P(3),PDPHI(IPI),1.0 280 CONTINUE C TEST 3 WRITE(20,*) WRITE(20,*) 'Test 3 - Dirichlet - {\phi}=X^2-Y^2+Z' WRITE(20,*) C Set up particular alpha and beta functions for Dirichlet/Neumann test C and type of boundary condition DO 340 ISP=1,NSP P(1)=SELCNT(ISP,1) P(2)=SELCNT(ISP,2) P(3)=SELCNT(ISP,3) SALPHA(ISP)=1.0D0 SBETA(ISP)=0.0D0 SF(ISP)=P(1)**2-P(2)**2+P(3) 340 CONTINUE DO 360 ISP=1,NSP SIPHI(ISP)=ZERO 360 CONTINUE DO 370 IPI=1,NPI PDIPHI(IPI)=ZERO 370 CONTINUE CALL LIBEM3(MAXNV,NV,VERTEX,MAXNS,NS,SELV, * MAXNPI,NPI,PINT, * SALPHA,SBETA,SF,SIPHI,PDIPHI, * LSOL,LVALID,EGEOM, * SPHI,SVEL,PDPHI, * L_SS,M_SSPLUS,L_DS,M_DS, * PERM,XORY,C,WKSPC1,WKSPC2,WKSPC3) WRITE(20,*) 'Solutions at the interior points' WRITE(20,*) ' P(1) P(2) P(3) Computed Exact' DO 380 IPI=1,NPI P(1)=PINT(IPI,1) P(2)=PINT(IPI,2) P(3)=PINT(IPI,3) WRITE(20,999) P(1),P(2),P(3),PDPHI(IPI),P(1)**2-P(2)**2+P(3) 380 CONTINUE 999 FORMAT(5F10.4) CLOSE(20) END REAL*8 FUNCTION PHI(P) REAL*8 P(3) PHI=P(1)*P(1)-2*P(2)*P(2)+P(3)*P(3) END REAL*8 FUNCTION DPHI(P,NP) REAL*8 P(3),NP(3) DPHI=NP(1)+NP(2)+NP(3) END